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ABSTRACT 

Aims. To model the interaction of the solar wind with the plasma tail of a comet by means of 
numerical simulations, taking into account the effects of viscous-like forces. 
Methods. A 2D hydrodynamical, two species, finite difference code has been developed for the 
solution of the time dependent continuity, momentum and energy conservation equations, as ap- 
plied to the problem at hand. 

Results. We compute the evolution of the plasma of cometary origin in the tail as well as the 
properties of the shocked solar wind plasma around it, as it transfers momentum on its passage 
by the tail. Velocity, density and temperature profiles across the tail are obtained. Several models 
with diff^erent flow parameters are considered in order to study the relative importance of viscous- 
like effects and the coupling between species on the flow dynamics. Assuming a Mach number 
equal to 2 for the incident solar wind as it flows past the comet's nucleus, the flow exhibits three 
transitions with location and properties depending on the Reynolds number for each species and 
on the ratio of the timescale for inter-species coupling to the crossing time of the free flowing so- 
lar wind. By comparing our results with the measurements taken in situ by the Giotto spacecraft 
during its flyby of comet Halley we constrain the flow parameters for both plasmas. 
Conclusions. In the context of our approximations, we find that our model is qualitatively con- 
sistent with the in situ measurements as long as the Reynolds number of the solar wind protons 
and of cometary H2O-1- ions is low, less than 100, suggesting that viscous-like momentum trans- 
port processes may play an important role in the interaction of the solar wind and the plasma 
environment of comets 

Key words, comets - solar wind — Hydrodynamics — comets: Halley - comet Giacobinni- 
Zinner - methods numerical 
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1. Introduction 



The nature of the interaction of the solar wind with the plasma environment of comets as they 
approach the Sun, has been under investigation since the early days of space physics as a discipline 
(Biermann 1951, Alfven 1957, see reviews by Cravens & Gombosi 2004, and Ip 2004). The basic 
elements of the interaction were developed in the 20 years following the work of Biermann (195 1), 
who proposed that the interaction between the solar wind and the comet's plasma is responsible for 
the observed aberration angle of plasma tails with respect to the Sun-comet radius vector. Based on 
the inefficiency of Coulomb collisional processes in the coupUng of the solar wind and cometary 
plasmas, Alfven (1957) proposed that the interplanetary magnetic field (IMF) is a fundamental 
ingredient in the solar wind-comet interaction; being responsible for channelling the cometary ions 
as it drapes into a magnetic tail. Biermann et al. (1967) suggested that as cometary ions are created 
and incorporated (picked-up) into the solar wind, the loading of the flow with this additional mass 
results in a modification of the flow properties as the solar wind approaches a comet; an idea further 
developed by Wallis (1973) (for a review see Szego et al. 2000). The IMF and mass loading are 
thus the main dynamical agents generally considered when developing models for the interaction 
of the solar wind with cometary ionospheres, as well as with other solar system bodies having an 
ionosphere and without a strong intrinsic magnetic field. 

However, as has been pointed out by Perez-de-Tejada et al. (1980) and Perez-de-Tejada (1989), 
several features of the flow dynamics in the cometosheath and plasma tail of comets can be at- 
tributed to the action of viscous-like forces as the solar wind interacts with cometary plasma. Such 
interaction processes are believed to be similar to those known to be occurring in other solar system 
bodies that have an ionosphere and no significant intrinsic magnetic field, particularly Venus and 
Mars (for a review see Perez-de-Tejada 1995, Perez-de-Tejada 2009 and references therein). In situ 
measurements indicate that, as in Venus and Mars, the solar wind flow in the ionosheath of comet 
Halley exhibits an intermediate transition, also called the "mystery transition", located approxi- 
mately half-way between the bow shock and the cometopause (Johnstone et al. 1986, Goldstein 
et al. 1986, Reme 1991, Perez-de-Tejada 1989 and references therein). Below this transition, as 
we approach the cometopause, the antisunward velocity of the shocked solar wind decreases in a 
manner consistent with a viscous boundary layer (Perez-de-Tejada 1989). Also indicative of the 
presence of viscous-Uke processes is that the temperature of the gas increases, and the density 
decreases, as we move from the intermediate transition to the cometopause. Taking the distance 
between the intermediate transition and the cometopause as the thickness of a viscous boundary 
layer, which depends on the eff^ective Reynolds number of the flow {Refs), Perez-de-Tejada (1989) 
estimated that R^s ~ 300 for the solar wind flow in the cometosheath is necessary to reproduce the 
flow properties measured in situ by the Giotto spacecraft on its flyby of comet Halley. 

An additional argument suggesting the importance of viscous-like effects in the dynamics of 
the flow in the cometosheath and tail regions, follows from the comparison of the magnitude of 
the terms corresponding to momentum transport due to viscous-like forces and J x B forces in 
the momentum conservation equation. Perez-de-Tejada (1999, 2000) has argued that downstream 
from the terminator in the ionosheath of Venus, a scenario analogous to the one considered in this 
paper, the fact that the flow is superalfvenic, as found from the in situ measurements of the Mariner 
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5 and Venera 10 spacecraft, suggests that viscous-like forces may dominate over J x B forces in 
the flow dynamics in the boundary layer formed in the interaction of solar wind and ionospheric 
plasma. If the flow is characterized by a low effective Reynolds number, R^s, this layer extends 
over a significant portion of the ionosheath of the planet. 

In comets, in situ measurements obtained during the passage of the ICE spacecraft through the 
tail of comet Giacobinni-Zinner (Bame et al. 1986, Slavin et al. 1986, Meyer- Vernet 1986, Reme 
1991) indicate that along the inbound trajectory (which lies slightly tailward of the comet nucleus) 
the magnetic field in the so-called transition and sheath regions, is approximately 10 nT, the number 
density is approximately 10 cm"-' and the tailward flow velocity varies from ~ 400 km/s (near the 
bow shock) down to 100 km/s. According to Perez-de-Tejada (1999), the ratio of viscous-like to 
magnetic forces is essentially the square of the alfvenic Mach number, = V^KB^/Snp). From 
the data of the ICE spacecraft cited above, we find that ranges between 4 and 40 across the 
cometosheath and hence, viscous-like stresses may dominate over J xB forces by a similar amount, 
or more, throughout the cometosheath region tailward of the nucleus. In the vicinity of the plasma 
tail, the measurements of the ICE spacecraft (Bame et al. 1986, Slavin et al. 1986) indicate that the 
midplane density, dominated by cometary ions, reaches values of 200 cm~^ at the point where the 
magnetic field is a maximum 50 nT. With flow speeds of approximately 20 km/s, the square of the 
alfvenic Mach number reaches a minimum value of 2-3 so that, even in the plasma tail, viscous-like 
forces are, at least, as important as J xB forces following the arguments of Perez-de-Tejada (1999). 

The fact that >> 1 in the cometosheath means that the magnetic energy density is much 
smaller than the kinetic energy associated with the inertia of the plasma. This implies that J x B 
forces are not the dominant dynamical factor responsible for the large scale properties of the flow 
in the region. In fact, one can argue that the formation of a magnetic tail is an indication that in the 
cometosheath, the large-scale magnetic field does not dominate the dynamics, it is merely carried 
around by the superalfvenic flow. If the dynamics were controlled by the magnetic forces, field fines 
would not bend onto a magnetic tail and the direction of the ion tail would not be essentially in 
the direction of the local solar wind velocity. We befieve that the magnetic field does play a crucial 
role in the momentum transfer between the solar wind and the cometary plasma, but it is the small 
scale, "turbulent" magnetic field component, that mediates the microscopic interaction between 
charged particles leading to the transfer of momentum that we are modelhng as an effectively 
viscous process. 

1.1. On the origin of "viscosity" 

The precise origin of the viscous-like momentum transfer processes invoked in the viscous flow 
interpretation of the intermediate transition, in the ionosheath of comet Halley and in other iono- 
spheric obstacles to the solar wind, is not yet clear. Typical properties of solar wind and come- 
tosheath plasma result in a "normal" viscosity, as it appears in the Navier-Stokes equations when 
derived from Boltzmaim's equation, that can be considered negligible in the flow dynamics. Using 
for example properties of the shocked solar wind in the vicinity of the tail measured at comet 
Giacobini-Zinner, n, = 10 cm^^ \B\ = 10 nT and T = 3 x 10^ K (Bame et al. 1986, Slavin et 
al. 1986) one calculates the viscosity coefficient resulting from particle interactions according to 
Spitzer (1962, eqn. 5-55) to be/i ~ 10'^ g cm"^ s ' . This extremely low value most likely repre- 
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sents a lower limit for the viscosity coefficient, since it reflects the abiUty to transport momentum 
across field fines in a plasma threaded by a strong, uniform magnetic field. A more appropriate 
expression for the plasma viscosity coefficient in the conditions of a cometosheath is probably 
given by the coefficient presented in Cravens et al. (1980), which corresponds to a plasma in a 
strongly fluctuating magnetic field. Perez-de-Tejada (2005) has calculated the viscosity coefficient 
for the ionosheath of Venus based on these results. If we use the same procedure to calculate the 
viscosity coefficient for the solar wind around the tail of a comet (with the conditions measured at 
Giacobini-Zinner) we find /i ~ 10 " g cm ' s"^ 

Wiffi typical values for the solar wind velocity and mass density in the cometosheath around 
the tail of comet Giacobini-Zinner, V = 200 km/s and p - 1.67 x 10 gm cm respectively 
(Bame et al. 1986), and adopting a characteristic lengthscale of 10^ km for the variation of the 
flow velocity (roughly the thickness of the sheath region), we find ffiat the corresponding Reynolds 
number for the flow, based on the "normal" viscosity coefficient estimated above, is Re > 10^ . This 
indicates that viscous effects resulting from the colfisions between particles in this environment 
are negUgible. Assuming that the Prandtl number is not very difl'erent from unity, as argued by 
Perez-de-Tejada (2005), we can also neglect heat conduction resulting from particle collisions. 

However, as in Venus and Mars, strong turbulence has been measured in the ionosheath of 
comets Halley and Giacobinni-Zinner (Baker et al. 1986, Scarf et al. 1986, Klimov et al. 1986, 
Tsurutani and Smith 1986) and, as it generally occurs in many fluid dynamics applications, turbu- 
lence is characterized (sometimes even defined) by a dramatic increase in the efficiency of transport 
processes, viscosity included, in the flow (Lesieur 1990). The Ukely importance of turbulent vis- 
cosity in this scenario is also expected in view of the large value of the Reynolds number estimated 
above. Also, as discussed by Shapiro et al. (1995) and Dobe et al. (1999 and references therein) 
conditions in the ionosheaffi of Venus and Mars favour ffie development of plasma instabifities lead- 
ing to efl'ective wave-particle interactions. If this mechanism operates also in the cometosheath, it 
may lead, as in these planets, to increased coupfing between the solar wind and cometary plasma 
in a viscous-Uke manner as suggested by Perez-de-Tejada (1989). In our opinion ffiis justifies a 
detailed study of the hypoffiesis of viscous-like effects on the flow dynamics in solar wind-comet 
interactions. It is the purpose of this paper to begin these investigations. 

In this paper we present results of 2D hydrodynamical, numerical simulations of the flow of 
solar wind and cometary H2O+ ions in the tail and tailward cometosheath of a comet. This is our 
first attempt to model the interaction of the solar wind with the plasma environment of a comet 
taking into account viscous-like forces which. We review the estimation of the effective Reynolds 
number of Perez-de-Tejada (1989), based on the comparison of in situ measurements at comet 
Halley with results from numerical simulations of the viscous-like, compressible flow of the solar 
wind over a dense, cold and slow velocity gas representing the plasma tail of a comet. We also study 
the relative importance of viscous-fike forces and the coupling between the fast moving protons of 
the solar wind and the slow H2O-1- ions in the tail. We do the latter by comparing models with 
different values of the effective Reynolds number, the parameter controlfing viscous-Uke effects, 
and the effective coupfing timescale between both species. 

The paper is organized as follows. In section 2 we present the formulation of the problem, the 
basic equations, approximations and parameters. Section 3 presents results of a series of simula- 
tions with different model parameters. A comparison of our results with in situ measurements at 
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comet Halley is discussed in section 4. 
present our conclusions. 
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Finally, in section 5 we summarize our main results and 



2. Formulation of the problem 

We model the interaction of the solar wind with the plasma tail of a comet using a 2D hydrody- 
namic, two species (a and b), finite difl'erence code that is an extension of the single species version 
presented in Reyes-Ruiz et al. (2008). Included in the dynamical equations is a coupUng term be- 
tween both species: solar wind protons and cometary ions, which we assume to be H2O-I- ions. This 
term allows the solar wind flow to get mass loaded with cometary ions as they diffuse upwards from 
the tail, and cometary ions to be accelerated by the fast, streaming solar wind. The coupling term 
is taken from the work of Szego et al (2000) who describe the treatment of mass loaded plasmas. 
However, in order to isolate the effects of the viscous-like forces, we do not consider the ongoing 
creation of new ions in the flow, by photoionization or any other mechanism, as is usually done in 
mass loading studies. Considering that we are modelhng only the flow in and around the tail of the 
comet, the only source of additional ions in our problem is through the boundary condition at the 
left hand edge of our simulation box (see §2.3). It is clear that the 2D character of our simulations 
is an approximation to the real problem and may not allow us to study some processes that may be 
essential for the dynamical evolution of the flow. We make this approximation considering that this 
is the first approach to the problem in which viscous-Uke forces are taken into account. We also 
neglect the effect of the IMF entrained in the solar wind flow, and leave for future work the study 
of the dynamical effects of J x B forces, although we do not expect these to be dominant in the 
region (see arguments in the Introduction section). 

Since we are interested in the gas dynamics in the tail we focus on the region behind the coma, 
starting from a few times 10'' km behind the comet's nucleus and extending downstream to a few 
times 10^ km as illustrated in Figure 1. 

2.1. Basic equations 

The present code solves the Euler equations for mass, momentum and energy conservation, includ- 
ing terms representing the viscous-like effects and interspecies coupUng due to turbulence and/or 
wave-particle interactions. In Cartesian coordinates and in conservative form, for species a, these 
can be written as: 




dU" dE" 
dt dx 



dy 



'ab 



(1) 
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Fig. 1. Illustration of the computational domain we use for our simulations. The box provides an 
approximate scale of the simulated region. Image of comet Halley taken the day of the encounter 
with Giotto by F.Miller, University of Michigan/CTIO (Brandt et al. 1992). 



and the inter-species coupUng term 

O 



(4) 



P^vativ; - V^) 



(5) 



In the preceding equations p" is the mass density of gas a, and V" are its velocity compo- 
nents, is its temperature and is the total energy density of species a defined by; 

It is important to point out that this form of the interspecies coupling term, although widely 
used in multispecies gas modelUng in various astrophysical scenarios (e.g. Schunk & Nagy, 1980, 
Draine, 1986, Cravens, 1991, Falle, 2003, Van Loo et al. 2009, Szego et al. 2000 and references 
therein), can be derived strictiy from the Boltzmann colUsion integral only for the case correspond- 
ing to Maxwell molecules (see for example Gombosi, 1994). We use it for lack of a similarly 
simple, alternative expression for charged particles, and must be considered an approximation of 
uncertain validity in our case. Schunk (1977) has discussed the modifications to these expressions 
for interspecies coupling for electrically charged molecules and in future contributions we shall 
explore the effect of such modifications. In the present calculations we have chosen this approach 
to modelling multispecies flow, which follows the dynamics of each species separately, instead 
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of an approach following a single fluid, composed of many different species, in order to clearly 
disentangle the widely different properties (p, V, T, etc) of solar wind and cometary ions. 

The coupling between species represented by the term 5"* in equation ([T) is taken from the 
work of Szego et al (2000), and has the form of the traditional coupling resulting from binary 
collisions. The term Vab contained in S"'', reflects the effective result of all processes able to transfer 
momentum and energy from one species to another. Note that in the adimensional form of the 
equations v^h is actually tgVah, that can be viewed as the ratio of the ffow crossing time, tg - 
L/Vo, to the inter-species coupling timescale, 1 / Va/,. In order to preserve the symmetry between the 
coupling terms for both species, guaranteed by the identity p^Vab - p'^Vba, we scale Vat as and 
Vha as p" with a single proportionality constant, v,,, which we take as uniform and constant. In our 
present code, v,, enters as a parameter that can be varied to compare the importance of inter-species 
coupling to viscous-like forces. 

= p''e" + ^p"(r'f (6) 

with e" being the internal energy per unit mass. In equations (O and (|4|i, the coefficients (i - 1 , 5) 
are the foflowing combinations of dimensionless numbers and the adiabatic index for the gas, y": 

^' = 

ki = (8) 
k", = (/ - 1), (9) 

k4 - \a ^ (10) 

y" 

k" = — , (11) 

ert 

where the Mach number (Mq), the Reynolds number {R^^ and Prandtl number (f r°) for the flow 
of gas fl, are defined respectively as: 

Mo = ^, (12) 

^ so 

PqVqL 

Pr" = Ll-P. (14) 

4 

Quantities with subindex o are those used for the normalization of the flow variables and param- 
eters, the reference sound speed is defined as Cso - '<Jy"Po/Po, c", is the specific heat at constant 
pressure for gas a, and L is the normalization for the spatial coordinates. For simplicity we have 
assumed that the flow parameters, ji and k, are uniform and that ji" - p", k" - /c", yU* = ju* and 

The terms T"^, , and T,", in equations ([3]) and (|4| represent the components of the viscous-like 
stress tensor given by: 

n'=^^_f^, (15) 
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dV" dV" 



and 

= + (17) 



2 dV} 4 dV^ 

3 dx 3 dy 



As is done in multiple fluid dynamics applications (Lesieur, 1990), we use the Boussinesq 
hypothesis in writing the Reynolds stress tensor, i.e. we adopt a "standard" form for the relation 
between the viscous-like stress tensor and the large scale flow velocity, using an effective viscos- 
ity coefficient that encapsulates turbulent viscosity as well as the possible effect of wave-particle 
interactions (Shapiro et al. 1995) or any other plasma instabilities leading to an increased coupling 
between ions in these coUisionless plasmas. 

Also, in equations (O and (|4]l, q" and q^, are the components of the effective heat flux vector for 
species a (under the Boussinesq hypothesis): 

dT" 

^--17' '''' 

and 

dT" 

« = --. (") 

Furthermore, we have assumed throughout this work that both gases are ideal so that: 
1 P" 



j" - 1 p« ' 



(20) 



with the equation of state, p" - p"RT" . We have assumed that both the solar wind plasma and the 
cometary plasma, in the tail region, are characterized by an adiabatic index, y" = = 5/3. In the 
section 4 of the paper we present some results for - 1.25, and discuss the eff'ects of changing 
this property of the cometary plasma. 

An analogous set of equations and definitions are written for species b, and both set of equa- 
tions, coupled by the source term S"* in equation ([1]), are solved simultaneously. 



2.2. Numerical code 

The set of equations described above is discretized in space using 2nd order finite difi^erences, 
and is advanced in time using an explicit, 2nd order MacCormack scheme (Anderson 1995). The 
implementation of the scheme is an extension of that described in Reyes-Ruiz et al. (2008), but 
now with the additional source term S in the equations of motion. In MacCormack's scheme the 
solution is advanced over one timestep by a sequence of intermediate steps, the predictor and 
corrector steps. In the predictor step an intermediate solution (t/*) is calculated from the values of 
the physical variables, U] j, at a given time, f, and position, (x,, yj), according to: 

U;, = l/b - [El.^j - E[j\ - C2 [n,,, - f ;J + At Sij (21) 

where ci - At/ Ax, ci - At I Ay and E', F' and 5' are evaluated with V according to (O, (|4| and (|5]). 
This predicted solution is then corrected to obtain the solution at the next time, f + At, using: 
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-[u[j + uIj-c,{eIj-eu;) 



(22) 



where E*, F* and S* are computed from U* using (O, (|4| and Q. Further details of the implemen- 
tation of MacCormack's scheme are given in Reyes-Ruiz et al. (2008). 

A final upgrade to our previous code is the ability to handle some types of non-uniform, carte- 
sian grids. For the simulations done in this work, the grid is defined by a series of (x,, yj) coordinates 
for which the spacing is arbitrary. In our simulations the x, points are geometrically distributed from 
Xmin to Xinax with nx elements. The yj points are equispaced at the initial location of the tail (from 
y = to y = 1 having 30 gridpoints) and geometrically distributed from y = 1 to y = ymax- In both 
series the common ratio is 1.02. The 2nd order approximation for the x-derivative of a function / 
at X, can be easily obtained from the Taylor series expansion of the function at x,_i and x,+i, and is 
given by: 



where Ax, = x;+i -x,. An analogous expression exists for the y-derivative. This grid allows a higher 
resolution in the vicinity of the region of strong interaction, while putting the y - y„ax boundary 
sufficiently far to avoid numerical artifacts in our results. 

2.3. Initial and boundary conditions 

The solution for the flow is evolved from the following initial conditions. A dense, cold, slow 
moving plasma representing the tail is located between y = and y - I. Both H2O-1- and H-H 
ions are present in the tail, but with protons much less abundant than H20h-. Between y = 1.5 and 
y - ymsLx, the gas has the properties of a shocked, hot, fast moving solar wind that contains both 
H2O-1- and H+ ions, with the number density of protons 50 times greater than H2O-1-. In all the 
calculations presented here, we have adopted a value Mg = 2 for the Mach number of the shocked 
solar wind incident on our computational domain. This assumption is made based on the results 
of Spreiter & Stahara (1980) who computed the the gas dynamics of the flow of the shocked solar 
wind in the ionosheath of Venus. Spreiter & Stahara (1980) found that the flow is characterized 
by M = 2, as the solar wind crosses the terminator of the planet (the line separating the day and 
night sides) and heads tailwards. In comets, we take the terminator to coincide approximately with 
the location of the nucleus. Between y = 1.0 and y = 1.5 there is transition region where the flow 
properties change smoothly in an exponential manner from those in the tail to those in the solar 
wind. The initial density of each species is taken to be: 





(23) 




Ptaii if y < 1 
0.32psw ify > 1-5 



(24) 
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We assume both species are moving initially with the same velocity: 



V/it = 0) 



I 



Vtaii if y < 1 
Vsw if3'> 1.5 



(25) 



V;.'\t = 0) 



0. 



(26) 



In normalized quantities Vsw - 1 and Vtaii - 0.01 . For the results shown here we use, in normal- 
ized variables, psw = 1 and ptaii = 400. The local temperature of both species is assumed the same 
inside the tail, with r°jj - T^^j = Ttaii and Jsw - lOOrtaii, with = 1 in normalized units. Outside 
the tail, fory > 1 .5, cometary ions are injected with a temperature an order of magnitude lower than 
the streaming solar wind protons, T"^ - IQT^^ - T^^. This choice of temperatures and densities is 
made to yield an initial pressure balance between the H20h- plasma (species b) inside the tail, and 
the proton plasma (solar wind, species a) outside. With p" - p"T" and = {m"plm!'p)p''T^ , m"^ and 
nip being the particle mass for species a and b respectively, we find that our choice of initial condi- 
tions is characterized by a pressure in-balance among each species. Whether the rapid movement 
of cometary ions resulting from this initial condition is prevented by the wrapped-around IMF over 
the comet's tail will be the subject of future studies. Although significantly different from the flow 
properties at later times in the simulations, for all the cases we have studied these initial conditions 
do not give rise to any long lasting instability in the flow so that the final state does not depend on 
their precise form or value. 

The boundary conditions are chosen to be consistent with the initial conditions. At the left 
boundary, x = Xmi„, the flow density and velocity follow exactly that given by the initial condition 
in equations (l23Tl-(l26b. Considering that the inflow to the comet's tail (y < 1) is subsonic, we 
allow the inflow pressure to float freely as a linear extrapolation of the active mesh values (e.g. 
Anderson 1995). The right side boundary, x - x^ax, corresponds to the commonly used outflow 
conditions for supersonic flows, namely the derivatives being zero for all flow variables. We have 
also run simulations with an outer boundary condition obtained from linearly extrapolating the flow 
variables, resulting only in minor diff'erences in the last gridpoints before the x = x^m- boundary. 



We have performed a series of simulations with different set of parameters R"'^ and Vo to determine 
the effect of viscous forces and inter-species coupling in the flow dynamics. For all cases consid- 
ered, the flow evolves from the prescribed initial condition [eqns. (I23l l-(l26]ll. passing through a fast 
transient phase, during which a considerable portion of the mass originally in the tail is eroded by 
the solar wind exiting our simulation domain. The relevance of this transient phase, lasting a few 
tens of solar wind crossing times (to = L/Vg), in relation to observed features in the evolution of 
the ion tail, will be analysed in a future publication. In this work we concentrate on the following, 
quiescent stage of evolution since, given its longer timescale for existence, is more likely to be 
encountered. In all cases, we present results for the flow velocity, density and temperature after a 
time long enough that a quasi-steady state has been reached. All results are presented in terms of 
normalized quantities as defined in the previous section. For a particular application, appropriate 
values of L, ¥„, Po and To can be chosen as exemplified in Section 4 for comet Halley. 



3. Results 
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Protons 




Fig. 2. Density contours (shades of gray) and flow geometry (velocity vectors) for Case 1 {Re"-'' = 
30, Vo =0.1) after 1234 simulation time units. The top panel shows the configuration for the proton 
plasma (species a) and the right side panel shows the "equiUbrium" configuration for cometary 
H2O-I- ions. Density and velocity are in normalized units. 

To determine the appropriate value of the effective Reynolds number for each species we con- 
sider the following. According to Perez-de-Tejada (1989) the geometry of the flow, measured in situ 
by the Giotto spacecraft in its fly-by comet Halley in march 1986, implies an effective Reynolds 
number around 350 for the shocked solar wind flow above the cometopause along the spacecraft 
trajectory. In contrast, in a similar region in the ionosheath of Venus, Perez-de-Tejada (1999) and 
Reyes-Ruiz et al (2008) estimate a value of the Reynolds number an order of magnitude smaUer 
(/?eff = 20), based on a comparison of in situ measurements (by the Venera 10 and Mariner 5 space- 
craft) at Venus with the flow properties derived from a numerical simulation of the viscous-hke 
solar wind-ionosphere interaction. To assess the estimation of Perez-de-Tejada (1989) we have 
conducted simulations with 3 different values of the Reynolds number. A high value, R^s = 100, 
similar to that estimated by Perez-de-Tejada (1989) for comet Halley; an intermediate value, /?eff 
= 30, comparable to the value estimated by Reyes-Ruiz et al. (2008) for the solar wind flow in the 
ionosheath of Venus; and a low value, R^s = 10, used to verify the tendency in the results as R^s is 
decreased. 

In most cases, the value of the effective Reynolds number for both species is assumed to be the 
same. In our view, the lack of knowledge of the precise mechanisms giving rise to the effective vis- 
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Fig. 3. Vertical profiles of the flow properties for Case 1 at three different positions; x - 2 (left 
column of panels), x = 5 (middle column) and x = 8 (right column). In all cases, gray lines 
indicate the properties of the proton plasma (species a) and black lines denote the properties of the 
H2O+ plasma (species h). The top row shows the x component of velocity, V"'*, the middle row 
shows the temperature, J"'^, and the bottom row shows the mass density, p"'*. AH quantities are in 
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cosity in these plasmas justifies this assumption. However, we have analysed a case with different 
values of the effective Reynolds number for each species in the Discussion section. 

The value of is also varied to explore the relative importance of inter-species coupling, 
versus viscous forces, which are proportional to R"^^. We will show results for 3 different cases: a 
strong coupling case characterized by Vo = 1, which can be interpreted as having the timescale for 
inter-species coupling equal to the solar wind crossing timescale, tg - LjY„\ a medium coupling 
case, in which the coupling timescale is an order of magnitude greater than the crossing timescale, 
Vo = 0.1; and a weak coupling case, for which Vo - 0.01, so that inter-species coupling effects are 
much smaller than other dynamical effects. 

To compare the state of the flow at the same time in its evolution for all cases, starting from 
the same initial condition, we have chosen, arbitrarily, to show results at f = 1234, with time units 
in multiples of the solar wind crossing time. The number of timesteps required to reach this time 
depends on the model parameters, for most cases less than 200 000 timesteps are required. 

3. 7. Effecf of inter-species coupling 

Our fiducial model. Case 1 , is characterized by model parameters = 30 and Vg = 0. 1 . In Figure 
|2] we show density contours and the flow velocity for each species. Initially the tail contained a 
uniform density, p" = 10 (protons) and = 400 (H2O+ ions) for y < I, and after t - 1234 ?„, 
a significant portion of the tail has been eroded by the effect of viscous forces and inter-species 
coupling. A shock wave is evident in the deflection of the flow velocity from the initial uniform 
distribution imposed by the boundary condition at x = Xj^m- Also noticeable is the strong velocity 
gradient around y = 2 which corresponds to the viscous boundary layer Both effects are also shown 
in Figure [3] where vertical profiles of the x component of velocity, V^, temperature, T, and mass 
density, p, are shown for three different x-positions, x = 2, 5 and 8; in the left, middle and right 
columns of each figure, respectively. 

In Figure [3] the shock front and the boundary layer can be identified at all 3 positions, but 
they are well separated only for x = 5 and x = 8, shown in the middle and right hand columns, 
respectively. For a given Vx profile, the shock front corresponds to the uppermost decrease from 
the uniform velocity {Vx = 1) in the free flowing solar wind. In the middle panel, corresponding to 
X = 5, this transition is located approximately at y = 5. A second transition, located approximately 
aty - 2.5 for x = 5, marks the top of the viscous boundary layer, below which the velocity drops 
sharply to the very low flow velocities in the middle of the tail. The shock front can also be seen 
as an increase in both temperature and density in the corresponding panels for each position. The 
temperature increase and density decrease characteristic of viscous boundary layers and found in 
previous studies of viscous flow over a flat plate (e.g. Reyes-Ruiz et al. 2008), is also observed in 
other cases modeled here. This clearly indicates that the region around y = 2 (at x = 5) is indeed a 
viscous boundary layer. 

For Case 2 we use the same Reynolds number, - 30, as in Case 1, but increase the impor- 
tance of inter-species coupling by using Vg - 1.0. A Figure showing the general flow geometry is 
not shown since no appreciable differences are found with Case 1 (shown in Figure|2]i. However, the 
vertical profiles of flow properties, shown in Figure|4] clearly illustrate the effect of a much stronger 
inter-species coupling used in this model. Namely, as both species are more tightly coupled, their 
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Fig. 4. Same as in Figure 3 but for simulation Case 2 (R^^ = 30, Vo = 1.0). 
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Fig. 5. Same as in Figure 3 but for simulation Case 3 (R^^ = 30, Vo = 0.01). 
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Fig. 6. Same as in Figure 2 but for simulation Case 4 (R"^g - 10, Vg - 0.1). 



velocity and temperature distribution tend to be almost identical. The density distribution conforms 
to the different boundary conditions for each species, since these are different, there is no reason 
why both densities should tend to equalize and they do not. Figure |4] shows that the shock front and 
the boundary layer are not well separated at the rightmost position shown, x - 2. From the shock 
front height and boundary layer thickness shown in the middle and right columns of Figure H] we 
see that both are proportional to the inter-species coupling (see §4). The shock front height at jc = 5, 
for example, changes from y = 4.7 for Case 1, to y = 5.5 in this case, while the thickness of the 
boundary layer goes from y - 2.8 to y = 3.2 as we increase the inter-species coupling parameter 
from 0.1 to 1 . 

In Figure|5]we show the results for Case 3 characterized by a very weak inter-species coupling, 
Vo = 0.01. The general flow geometry (not shown) is very similar to that in Figure|2] A comparison 
of Figure |5] (weak coupling) with Figures |3] and |4] (medium and strong coupling respectively), 
clearly shows that in Case 3 the dynamics of both species is essentially uncoupled. The location 
of the shock front and the top of the boundary layer are different for each species. For example, at 
X - 5, only for the cometary H20h- ions the shock front and boundary layer are clearly separated. 
For the H2O-1- ions the shock front is located approximately at y = 5 and the top of the boundary 
layer is at y - 2.5, while for protons the shock front and the top of the boundary layer are both 
located around y = 2. 
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To analyse the effect of the viscous-like momentum transport between the solar wind and material 
in the comet's plasma tail, we compare 3 simulations with the same inter-species coupling param- 
eter, Vo = 0.1, but different values of the effective Reynolds number. Figures |6] and [T] shown the 
resulting flow geometry and vertical profiles, respectively, for our Case 4, characterized by a higher 
viscosity corresponding to a lower effective Reynolds number, R^ft -10, than Case 1. Comparing 
the global geometry of the flow in this case (Figure |6]) with that in a case with greater Reynolds 
number, R"'^ - 30 (Figure|2]i we see that after 1234 crossing times, the erosion of the tail is much 
greater in this high viscosity case for both species. This result is expected as well as the increase in 
the thickness of the boundary layer as we decrease the effective Reynolds number. This is clearly 
seen when comparing the vertical profiles of the flow properties shown in Figure [3] (medium vis- 
cosity) and Figure |2] (high viscosity). For example, as shown in Figure |7] for x - 5, the top of the 
boundary layer increases from y = 2.8 for Re"''' = 30 to approximately y = 3.7 for R"J = 10. 
The increased thickness of the boundary layer as we decrease R"J^, effectively represents a more 
blunt obstacle to the solar wind flow. Hence, the height of the boundary layer also increases as we 
decrease Re"-''. This is also shown in Figure Q where, for example at jc = 5, the height of the shock 
front is located approximately at y = 7; about 2 scale units higher than the shock front location for 
the model with lower effective viscosity (Figure[3]l. 

The tendency seen in going from high (R''J^ - 10) to medium effective viscosity (R"J^ - 30) is 
confirmed by comparing with results with an even smaller viscosity, such as Case 5 which corre- 
sponds to a model with Re"''' - 100, shown in Figures[8]and|9] As expected, a decreased viscosity 
leads to significantly less erosion of the tail than in Cases 1 and 4 (medium and high viscosity 
respectively) as shown in Figure [8] Also, as discussed above and as shown in Figure |9] the top of 
the boundary layer decreases as we increase the Reynolds number, and consequently the location 
of the shock front also decreases. For example, in the profiles corresponding to jc = 5 in Figure 
|9l we find that the top of the boundary layer decreases from y - 2.8 for R"'g - 30 to y = 2.2 for 
^eff ~ regards to the location of the shock front, this goes from y = 4.8 for R"J = 30 to 

approximately y = 3.7 for^^j^* = 100. 

4. Discussion 

In view of the uncertainty about the precise physical mechanisms giving rise to the effective vis- 
cosity, we have assumed that the effective Reynolds number for both species is the same in the 
calculations presented above. However, we have also carried out simulations having distinct ef- 
fective Reynolds number for each species and find that for a medium value of the inter-species 
coupling, Vo =0.1, the results are almost identical to those with a single value for the effective 
Reynolds number for both species (equal to the effective Reynolds number of species b). For ex- 
ample, for a case with 7?"^^ = 100, R'^g =10 and v,, = 0.1, the vertical profile of flow properties for 
the H2O-1- ions at all x locations is almost identical to that shown in Figure|2]for Case 4, character- 
ized by R"g - R'^g =10 and = 0.1. The vertical profile of flow properties for the protons, while 
not identical, is still very similar to Case 4. This suggests that viscous stresses, particularly in the 
species that dominates the mass of the problem, are the dominant factor in the flow dynamics. 
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Fig. 8. Same as in Figure 2 but for simulation Case 5 (Rl^ - 100, Vo - 0.1). 

As mentioned in section 2, in the results presented above we have assumed that the adiabatic 
index for both species is the same, y" - y'' - 1 .67. While this value of y can be safely assumed for 
the solar wind plasma (assuming thermal equilibrium for the species), it is not so clearly valid for 
the H20^ plasma in which the excitation of rotational and vibrational degrees of freedom may lead 
to a lower value of y (again assuming thermal equilibrium for the species). In order to illustrate 
the effect of a different, lower value of the adiabatic index for cometary plasma, we have also con- 
ducted simulations with a value y'' - 1 .25 for the adiabatic index of the H20^ plasma. This value 
corresponds to a gas composed of triatomic molecules in thermal equilibrium at a high enough 
temperature for all molecular degrees of freedom to be excited. Results for this case, y" - 1 .67 
and - 1.25, with the same effective viscosity and interspecies coupling parameters as Case 1 
(R"g - 30 and Vg = 0.1) are shown in Figure [10] which shows the vertical profiles of V^-, T and p 
for both species in both cases. 

Clearly evident when comparing Figure [10] (y'' - 1 .25) and Figure [3] (y'' - 1 .67) is the fact 
that if the cometary plasma is characterized by a lower value of the adiabatic index, the heating of 
the HaO^ plasma is significantly reduced in the boundary layer, since part of the dissipated energy 
goes to the excitation of the additional degrees of freedom corresponding to the lower value of y. 
This leads to less plasma expansion in the region and a thinner velocity boundary layer. The height 
of the shock front is consequently reduced. In future contributions we shall address the issue of the 
appropriate value of y for the cometary plasma. 
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Fig. 9. Same as in Figure 3 but for simulation Case 5 (R^^ = 100, Vo = 0.1). 
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Fig. 10. Comparison of the vertical profiles of flow properties for a case characterized by the same 
value of /?etf and Vo as Case 1, but with - 1.25. Profiles at x = 2 (left column of panels), 
X = 5 (middle column) and x = 8 (right column) are shown. Gray lines indicate the properties of 
the proton plasma and black lines denote the properties of the H2O+ plasma. All quantities are in 
normalized units. 
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4.1. Comparison with in situ measurements 



A comparison of our results with the in situ measurements made by the Giotto spacecraft, as it 
flew by comet HaUey in March of 1986, is not straightforward. The simplified geometry we are 
using in our simulations to study the interaction in the tail region exclusively, precludes a direct 
comparison. Nevertheless, some insight into the implications of our results can be obtained from a 
simplified comparison. 

Once a particular application scenario has been chosen, values for the characteristic length, ve- 
locity, density and temperature used in the adimensionalization of the equations of motion (section 
|2l) can be established. For comet Halley, using the in situ measurements reported in Goldstein et al. 
(1986), Johnstone et al. (1986) and Perez-de-Tejada (1989), we adopt L = 150,000 km, V„ = 250 
km/s, p„ = 1.67 X 10-23 gm/cm^ and T,, = 2.5 xlO^ K. 

According to Johnstone et al. (1986), the Giotto spacecraft observed 3 distinct transitions in 
the plasma properties on its inbound trajectory towards comet Halley 's nuclear region: (1) The 
outermost transition occurs about 900000 km from the point of closest approach and can be iden- 
tified as the bow shock crossing. (2) The cometopause, where the density of cometary ions sharply 
increases, can be located at around 150000 km from closest approach (Perez-de-Tejada, 1989). 
(3) Approximately midway between the shock location and the cometopause, at about 400000 km 
from closest approach, the so called intermediate transition signals the top of the viscous bound- 
ary layer according to the viscous flow interpretation of the solar-wind-comet interaction given by 
Perez-de-Tejada (1989). Pending a more detailed comparison of the Giotto measurements with the 
results of our simulations, which should take into account the full geometry of the problem, let us 
identify the cometopause detected in the measurements with the region of very strong H2O4- den- 
sity increase in our simulations, located aty - 1.0 (approximately) in our normalized units. Under 
this assumption, in Figure [TT] we compare the thickness of the boundary layer and the height of 
the shock front evaluated from our simulation results at x = 5, for models with different effective 
Reynolds number {R"^) and inter-species coupling parameter (v^). As akeady seen, both the thick- 
ness of the boundary layer and the height of the shock front decrease with increasing Reynolds 
number so that, almost irrespective of the value of v„, a low value of R"^ is required to explain the 
measured transition locations. 

Also evident in Figure [TT] is the dependence of the transition locations on the value of the 
inter-species coupling parameter In simulations with a strong inter-species coupling, the solar 
wind ions are able to transfer momentum to cometary ions more efficiently giving rise to a thicker 
boundary layer and higher shock front. The opposite is true when both species are weakly coupled 
(Vo = 0.01). In such case solar wind ions flow by cometary plasma interacting very weakly. Less 
momentum is transferred between the solar wind and cometary plasma in a situation reminiscent 
of a high Reynolds number case. Our analysis of scale-heights is based not only on the properties 
of the velocity profiles in our simulations. As pointed out by Perez-de-Tejada (1989), there are 
corresponding changes in the density and temperature of the gas as one enters a boundary layer 
The heating and expansion characteristic of viscous boundary layers are also found in our results, 
particularly for cases with low Reynolds number and high inter-species coupling parameter. 

It is worth mentioning that somewhat similar properties of the flow were also measured by the 
ICE spacecraft in its flyby through comet Giacobinni-Zinner as discussed in Ip (2004). A com- 
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parison of the location of the transition from the sheath region to the so-called transition region 
and the bow shock location as estimated by Reme (1991), corrected for the different height of the 
cometopause, yields very similar relative positions to those shown by the dotted lines in Figure 
[TT]for the transitions in comet Halley. In future work we will address the differences in the flow 
properties measured in comet Halley and in comet Giacobinni-Zinner. 

4.2. Implications for 3D geometry 

It is important to emphasize that the geometry presented in this paper is derived from a 2D model. 
In Venus, as discussed by Perez-de-Tejada (1995), the viscous-like interaction between solar wind 
and ionospheric plasmas takes place preferentially over the magnetic poles of the planet (defined in 
terms of the incident IMF), where the pile-up of magnetic field lines is less than around equatorial 
latitudes. According to Perez-de-Tejada (1995), up to about 80° SZA, the piled-up magnetic field 
over the dayside ionosphere and along the flanks, inhibits in some degree a direct, viscous-like 
interaction between solar wind and ionospheric plasmas. 

If we apply these ideas to the solar wind-comet interaction, this implies that the flow prop- 
erties we have computed here, correspond more closely to locations over and downstream from 
the magnetic poles of the comet. For different locations along the tail, the piled-up magnetic field 
may prevent an efficient viscous-like dragging of ionospheric material, J x B forces may be more 
important and the flow dynamics may be better modeled in terms of an MHD model as those of 
Wegmann (2002) and Jia et al. (2007). As the IMF is constantly changing direction on a wide 
range of amplitudes and timescales, the region of viscous-like interaction between the solar wind 
and cometary plasma, changes with time. Given the typical IMF orientation is approximately in 
the ecliptic plane, one should expect that the flow within +/- 20°, measured in the y-direction (as 
typically defined) from the magnetic poles of the comet, is best described by our model. 

5. Conclusions 

We have presented results for the numerical simulation of the interaction between the solar wind 
and the plasma in the tail of a comet, taking into account the effect of viscous-like stresses previ- 
ously argued to be important by Perez-de-Tejada et al (1980). To our knowledge, this is the first 
time that viscous-like effects have been incorporated into such studies. Our results indicate the ex- 
istence of 3 distinct transitions in the flow properties: outermost we find a shock front, innermost 
we have the cometopause and an intermediate transition which we can identify with the height of 
the boundary layer characterized by a fast decline in the anti-sunward flow velocity, and the onset 
of plasma heating and expansion due to viscous-like dissipation. The location of these transitions 
depends on the flow parameters, namely the effective Reynolds number of the flow for each species, 
R^g, and the inter-species coupling parameter, v„. 

By comparing the flow properties from our numerical simulations to the location of the shock 
front and intermediate transition, as measured by the Giotto spacecraft as it approached the nucleus 
of comet Halley, we find that, almost irrespective of the strength of the inter-species coupling, v,,; a 
low value of the effective Reynolds number, approximately R"^ < 20 for both species, is required 
to reproduce the measured transition locations. This implies, in the context of our model, that the 
measured flow properties cannot be explained if one does not take into account the viscous-like 
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forces in the interaction of the solar wind and the plasma tail of a comet. Although the conclusions 
drawn from this study are strictly appUcable only to comet Halley and solar wind conditions at the 
time the in situ measurements were taken, one may speculate that viscous-like processes may be 
important in the solar wind-comet interaction in general. 

It is important to emphasize that, this being the first attempt to include viscous-like forces in 
the numerical simulation of the interaction of the solar wind with a comet's plasma environment, 
there are many pending issues still to be addressed that could have potentially important conse- 
quences on the details of the solutions obtained under our simplified treatment. First and foremost, 
the precise forms we are using for the viscous like stress and efl'ective interspecies couphng, may 
be questioned. As we have argued in the Introduction, plasma properties imply that "normal" vis- 
cosity is negligible in the region under consideration. Hence, we are invoking an effective viscos- 
ity presumably resulting from plasma turbulence and/or wave-particle interactions. However, the 
precise form of the terms corresponding to viscous-Uke momentum transfer in the equations of 
motion (Bousinessq hypothesis) is not formally demonstrated. Also, as we have discussed in the 
Formulation section of the paper, the interspecies coupling terms we are using can not be strictly 
derived for a plasma as the one we are modelling. In view of these arguments, one may consider 
that the work reported in this paper is only an academic exercise of questionable applicability to 
the problem of solar wind-cometary plasma interaction. In such case, a similar conclusion must be 
reached in regards to many other studies of fluid dynamics that use similar approaches to modelling 
eff'ective viscosity and interspecies coupling. 

Additional important effects still to be considered are the following: geometrical effects due to 
the curvature of the ionosphere are required for a more direct, quantitative comparison between in 
situ measurements by the Giotto spacecraft and the results of simulations; the interaction of the 
charged species with neutral gas ejected from the comet which, especially in the vicinity of the 
nucleus, is the most abundant species; the effect of the magnetic field on the flow (particularly 
in the dayside and around the midplane of the near-tail region), 3D effects, incoming flow time 
dependence, etc. We believe that the further assessment of the relevance of these factors is beyond 
the present study. They are the subject of work currentiy in progress and wiU be reported in future 
contributions. 
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Fig. 11. Height of the shock front, i/j/, (gray lines with squares) and thickness of the boundary 
layer, 5 (black lines with triangles), as a function of effective Reynolds number, for a set of models 
with different value for the inter-species coupling parameter, Vg. Values for these scale-heights 
correspond to x = 5 in our model. The dotted lines indicate the height of the shock front (gray) 
and the location of the intermediate transition (black) during the inbound portion of Giotto's flyby 
through the tail of comet Halley. 



